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We consider RKKY interaction between two magnetic impurities in graphene. The considera- 
tion is based on the perturbation theory for the thermodynamic potential in the imaginary time 
representation. We analyze the symmetry of the RKKY interaction on the bipartite lattice at half 
fiUing. Our analytical calculation of the interaction is based on direct evaluation of real space spin 
susceptibility. 
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INTRODUCTION 

Since graphene was first isolated experimentally 
it is in the focus of attention of both theorists and ex- 
perimentalists. Many physical phenomena, well stud- 
ied in "traditional" solid state physics look quite dif- 
ferent in graphene. In this paper we will talk about 
the Ruderman-Kittel-Kasuya-Yosida (RKKY) interac- 
tion, first studied (in a normal metal) more than 60 years 
ago 0-0]. This interaction is the effective exchange be- 
tween two magnetic impurities in a non-magnetic host, 
obtained as the second order perturbation with respect 
to exchange interaction between the magnetic impurity 
and the itinerant electrons of the host. 

Quite a few theoretical papers published recently con- 
sidered RKKY interaction in graphene Though 
analysis of the RKKY interaction is simple in princi- 
ple, calculation of the integrals defining the interaction 
(whether analytical or numerical) can pose some prob- 
lems. However, substantial progress was achieved in the 
field. 

Our interest in the RKKY interaction in graphene 
started from learning about the theorem stating that for 
any half-filled bipartite lattice the exchange interaction 
between the magnetic adatoms is ferromagnetic, if the 
adatoms belong to the same sublattice, and antiferro- 
magnetic, if the adatoms belong to different sublattices 
Q. Also, in this paper and in the following one [lo| . 
in the approximation of the linear dispersion law for the 
electrons, the RKKY interaction in graphene was cal- 
culated analytically. However, the integrals obtained in 
both papers turned out to be divergent, and the compli- 
cated (and to some extent arbitrary) cut-off procedure 
was implemented to obtain from these integrals the finite 
results. So we started to look for the procedure which will 
allow to eliminate this problem. The second reason for 
our interest was the fact that the theorem, mentioned 
above, was challenged The claim was that the proof 
is based on calculation of the magnetic susceptibility of 
the free electron gas by the imaginary-time method, de- 
manding later analytic continuation from the imaginary 
frequencies to the real ones. On the other hand, consid- 
eration by the real-time method, presented in Ref. 12 1 



does not support the statement of the theorem. To clar- 
ify the situation and get rid of the shortages mentioned 
above, we decided to analyze the problem of RKKY in- 
teraction in graphene from the scratch. 



RKKY INTERACTION 

We consider two magnetic impurities at the sites i and 
j and assume a contact exchange interaction between the 
electrons and the magnetic impurities. Thus the total 
Hamiltonian of the system is 

Ht = H + H„u = H- JS,-s, - JSj-Sj, (1) 

where H is the Hamiltonian of the electron system. Si is 
the spins of the impurity and is the spin of itinerant 
electrons at site i. 

Our consideration is based on the perturbation theory 
for the thermodynamic potential [l3|. The correction to 
the thermodynamic potential due to interaction is 

An = -Tin (S) = -Tlntr js" • e-"^^/Z^ , (2) 

where the iS-matrix is given by the equation 

fl/T ~) 



S = exp ■ 



Hint{T)dT 



(3) 



Writing down in the second quantization representa- 
tion 

S« = ^cfaO-a/jCi/S, (4) 

the second order term of the expansion with respect to 
the interaction is 



(5) 



1/T pI/T 



^0 



Notice that we have ignored the terms proportional to Sf 
and S|, because they are irrelevant for our calculation of 
the effective interaction between the adatoms spins. 



2 



Leaving aside the question about the spin structure of 
the two-particle Green's function standing in the r.h.s. of 
Eq. ^ (for interacting electrons), further on we assume 
that the electrons are non-interacting. This will allow us 
to use Wick theorem and present the correlator from Eq. 
([5]) in the form 



Gl3fiiJ;Ti - T2)Qsa{i,i\T2 - Ti), (6) 



where 



QMi.hi,Ti - ra) = - (Tr |c,^(Ti)cj^(T2)|^ (7) 

is the Matsubara Green's function [l^. We can connect 
Qp^ with the Green's function of spinless electron 

Gl3jiiJ,Ti - Ta) = -Sp^ (Tr |c,(ti)c](t2)|^ . (8) 

Presence of delta-symbols allows to perform summation 
with respect to spin indices in Eq. (0) 



which gives 



(9) 



(10) 



where 



1 /-i/^ 

Xtj=--; G{hj;r)g{j,i;-T)dT (11) 



is the free electrons static real space spin susceptibility. 
Thus we obtain 



3-1 



(12) 



Eq. (Ilip was applied to calculation of RKKY interaction 
in graphene for the first time, to the best of our knowl- 
edge, in Ref. [lij . 

The Green's function can be easily written down us- 
ing representation of eigenvectors and eigenvalues of the 
operator H 



It is 



{H - En) Un = 0. 



(l-nFi^n)), T>0 



(13) 



r < ' 



(14) 



where £,„ = En — jJ, and nF{^) = (e^^ -f- 1) ^ is the Fermi 
distribution function. 



SYMMETRY OF THE RKKY INTERACTION ON 
THE HALF-FILLED BIPARTITE LATTICE 

In this Section we'll consider the Hamiltonian of the 
free electrons in tight-binding representation 



H — tijc\cj. 



(15) 



Bipartite lattice we'll understand in the sense, that all the 
sites can be divided in two sublattices, and there is only 
inter-sublattice hopping (no intra-sublattice hopping). 
Thus the Hamiltonian H in matrix representation is 



H 



T 

Tt 



(16) 



where T is some matrix NxM ( the first N sites belong to 
the sublattice A and the last M sites belong to sublattice 
B). 

Consider a matrix of even more general form than ([T| 



H 



OnxN BnxM 
CmxN Omxm 



(17) 



B and C are some arbitrary matrices. The spectrum of 
the matrix H can be found from a secular equation 



-EInxn 

CmxN 



BnxM 

-EImxm 



0. 



(18) 



In Ref. 15| it is proved the following property of the 
determinant of the block matrix 



CmxN DmxM 



= \A\ \D-CA-^B\ 



(19) 



which is valid, provided |^| ^ 0. For non-zero eigen- 
values of the matrix H , we can apply Eq. (|19p to the 
determinant to get 



\E^I. 



MxM 



-CB\ =0. 



(20) 



Thus the spectrum of the bipartite Hamiltonian is sym- 
metric, that is non-zero eigenvalues of the matrix H are 
present in pares [E, —E). 

If we write down Eq. (|13p explicitly in a matrix form 



-Enl T 

Tt -E^I 
it becomes obvious that 

Un{i) = ±U„(«), 



0, 



(21) 



(22) 



where u„ is the eigenfunction corresponding to En and 
Un is the eigenfunction corresponding to —En, and in the 
r.h.s. of Eq. ([^ there is plus sign if the site i belongs to 
one sublattice, and there is minus sign if the site belongs 
to the opposite sublattice. 

Eq. (l^^ and the fact that for /i = we have np{£^m) = 
1 — ni?(^„j), immediately convince us that the terms with 
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non-zero energy in Eq. ([T^ are pairwise antisymmetric 
(with respect to simultaneous transformation r — >■ — r, 
j J and complex conjugation) for the sites i and j be- 
longing to the same sublattice, and pairwise symmetric 
for the sites i and j belonging to opposite sublattices. 
The term (terms) with E = is antisymmetric with re- 
spect to the above mentioned transformation, no matter 
which sublattices the sites belong to. Thus for the sites 
i and j belonging to the same sublattice 



Gij,i;-T) = -g*{i,j;T). 



(23) 



For the sites i and j belonging to different sublattices 



(24) 



provided there are no zero energy states, or we can ne- 
glect there contribution to the Green's function. 

Thus for the case considered, Eq. (fTT|) gives ferro- 
magnetic exchange between magnetic impurities on the 
same sublattice and antiferromagnetic exchange between 
impurities on opposite sublattices (under the restriction 
presented above). 



ANALYTIC CALCULATION OF THE RKKY 
INTERACTION IN GRAPHENE 

In calculations of the RKKY interaction in graphene 
the ^„ in Eq. (fH|l turns into J d^P, where a is the 
carbon-carbon distance. (Actually, there should appear 
a numerical multiplier, connecting the area of the ele- 
mentary cell with a^, but we decided to discard it, which 
is equivalent to some numerical renormalization of J.) 
Also 



Un{i) 



(25) 



where ipp is the appropriate component of spinor elec- 
tron wave-function (depending upon which sublattice the 
magnetic adatom belongs to) in momentum representa- 
tion. 

Further on the integration with respect to cPp we'll 
treat as the integration in the vicinity of two Dirac points 
K, K' and present p = K(K') -f k. The wave function 
for the momentum around Dirac points K and K' has 
respectively the form 



V'i^,K(k) = 
V'i^.K' (k) 



72 V ^'^'^'^''^ 
1 / 



(26) 



where v = ±1 corresponds to electron and hole band [iq : 
the upper line of the spinor refers to the sublattice A and 
the lower line refers to the sublattice B. 

We'll assume that T = and the Fermi energy is at 
the Dirac points; E+Ck) and ii'-(k) would be electron 



and hole energy. Then Eq. (fTi|) takes the form: for i and 
j belonging to the same sublattice 



C/^^(z,j;r>0) = - 



1 r 



(2n) 



^2j^gikR.j-£;+(k) 



(27) 



and for i and j belonging to different sublattices 



c;^^(z,j;t>0) = - 



1 a' 



2 (27r)2 



,i(K+k)-Ri3-ifl 



(28) 



For T < we should change the sign of the Green's func- 
tions and substitute E- for E+. 

For isotropic dispersion law -B(k) — E{k) we can per- 
form the angle integration in Eqs. P7|) and (pSj) to get 

^2i^gZk.R.,-£;(fe)r ^ / dfcfcjp(/ji?)e-B(fc)^ (29) 
"'0 

(Jo and Ji are the Bessel function of zero and first order 
respectively, and 0r is the angle between the vectors K — 
K' and Ki/). 

For the linear dispersion law 



E±{k) = ±VFk, 
using mathematical identity jl7l | 



(30) 



(31) 



(-1) 



dp 



we can explicitly perform the remaining integration. Cal- 
culating integrals we obtain 



3a^ 



[l-Hcos((K-K')-Rij)] (32) 
1 - cos((K - K')-Ry - 29r)] . 



(33) 



The approach presented above can be easily applied 
to the bilayer graphene. We'll consider Bernal {A — B) 
stacking. Because the low-energy modes are localized 
on A and B sites 18|, we consider RKKY interaction of 
the magnetic adatoms siting on top of carbon atom in A 
and/or B sites. The low-energy modes are characterized 
by the spectrum 



i?±(k)=± 



2m 



(34) 
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and wave functions 



V2\ 



(35) 



where this time the upper hue of the spinor refers to the 
sublattice A and the lower Hne refers to the sublattice B 
[3] (we ignore the trigonal warping). So for the case of 
bilayer we reproduce Eq. (P7| (of course, the result for 
equally refers to G^^); Eq. ^ is changed to 



gi(K+k)-Ri3-2jefc _ gj(K'+k)-Rij+2iefc 



Calculation of would demand the integral [lE 



(36) 



Jo(x) exp(—px )xdx = — exp 
2p 



1 

Ap 



(37) 



After simple calculus we obtain for bilayer graphene 



,AA 



(R», 



167r2i?2 



[l + cos((K-K')-Ry)]. (38) 



We return to monolayer graphene. The case of mag- 
netic adatom siting on top of carbon atom, certainly does 
not exhaust all the pos sibilities for the adatom positions 
in graphene lattice [111 . [l4| . However, under rather gen- 
eral assumptions the specific position of the adatom can 
be taken into account by changing in Eq. (|14p the prod- 
uct of the components of the spinor wave function ijjp to 
an appropriate matrix element. Thus, using the results of 
Ref. [11], for the case of substitutional impurities instead 
of Eqs. (HH) and (HH]) we obtain 



Cy^^(^,J;r>0) = -- 



^AB 



iK Ri 



iK'Ri 



Cy-^^(z,j;r>0) = - 



1 



2 (27r) 



g4(K+k)-R,j-3iefc _ gi(K'-k)-R,j+3iefc 



(39) 



After simple calculus we obtain 

X^-'- (R„) = [1 + cos((K - K')-R„)] 



and 

(Ry) 



vfR'' 

XSaS_ 



X 



[l-cos((K-K')-Ry -60r) 



where X^-^^^ and X^^^^ can be easily calculated ana- 
lytically. 



DISCUSSION 

In this Section we would like to compare our results 
with the previously obtained ones and additionally justify 
our line of reasoning. 

The correction to the thermodynamic potential can be 
also written down using frequency representation [l3| . 
which would give 

n.m 



1 



1 



(40) 



where uj — ttT(21 + 1) {I is an integer) is Matsubara fre- 
quency. (Eq. p4p was taken into account.) Performing 
in Eq. (|4U|) summation with respect to Matsubara fre- 
quencies we obtain 



<(i)u„(i)u„(j)<„(j)- 



(41) 



This is in fact the result obtained originally 043 j by 
using standard quantum mechanics (off-the-energy shell) 
perturbation theory. 

We used the term bipartite lattice, but actually neither 
the symmetry of spectrum , presented after Eq. (|20p , nor 
the symmetry of wave functions presented in Eq. (|22p do 
not require any space periodicity (or any order at all) in 
the position of the sites. They even do not require that 
the Hamiltonian will be Hermitian, so they remain, say, 
in non-Hermitian quantum mechanics. 

The results of Ref. 1^ correspond to Eq. (j¥T|) with a 
small but substantial difference: the terms with ^„ — 
are discarded [1^, which breaks the symmetry of the 
RKKY interaction we discussed. We want now to con- 
sider a simple toy model to additionally explain that 
these diagonal terms are relevant and should be where 
they are. Our arguments will follow the consideration of 
the magnetism of electron gas in Ref. 



21| 



Let the spectrum of H consists of pairs of states having 
the same energy, and Hint has non-zero matrix elements 
only between the states belonging to the same pair. Then 
the quantum mechanical problem of finding the spectrum 
of the Hamiltonian Ht can be solved exactly, each dou- 

(12) 

blet is split, £"„ ' ' — i?n± |^n,i;n,2|- The thermodynamic 
potential is 



[En ± IK, 



(42) 



where Vf'{E) is the thermodynamic potential of the iso- 
lated level with the energy E. Expanding with respect 
to interaction we obtain 



2o0 



dEl 



n,l;n,2 — 



/ ^ Kl=, \Vn,l;n,2\ 



dEn 



(43) 
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which corresponds to the diagonal terms in Eq. (HT]) . The 
issue of diagonal terms can be also connected to the dif- 
ference between the real- and imaginary-time approaches 
the authors of Ref. emphasize in their paper. Our 
opinion is that calculation of magnetic susceptibility us- 
ing real-time method (Kubo formula) gives the adiabatic 
susceptibility. On the other hand, for the calculation of 
the RKKY interaction we need the isothermal suscepti- 
bility, which is given by the imaginary-time method. 

Analytical calculations of the RKKY interaction can 
be done using Eq. ([1T|) . In this case instead of Eqs. 
and (1331) we would obtain 



[l-fcos((K - K')-Rj. 



1 



dxxJo{x) I dx'x'Jo{x') — ^ (44) 
Jo Jq X + X 

(R„) = -T^^^ [1 - cos((K - K')-Ry - 20^)] 



1 



dxxJi{x) I dx'x' Ji{x')- ^. 
Jo X + X 



(45) 



Eqs. ([33]) and (|45p are particularly convenient to be com- 
pared with the results of Ref. [l^l ■ Using the identity (l7j 



X + z 



Jv{cx)dx 



2 cos VTT 



[H_„(cz)-y_,(cz)], 



(46) 



where H^(z) is the Struve function and Y^{z) is the Neu- 
mann function, we can present integrals in Eqs. (j44p and 
(133) as 



dxx^Joix) 



Fo(a;)-Ho(x) 



^ / dxx^Ji{x)[Y^i{x) - H_i(a; 
2 Jo 



)]■ (47) 



These integrals are similar to those standing in Eqs. (18) 
and (25) of Ref. 0, but contrary to the latter, our 
integrals diverge. This is guaranteed by the asymptotics 
of Struve functions 



H,(x)-y,(a;) 



X\ 

2. 



o{{x/2r 



(48) 



A deficiency of the previous analytic calculations of the 
RKKY interaction in graphene is, to our mind, not due to 
them using the frequency representation of the Green's 
function (though we find the imaginary time represen- 
tation more convenient for the calculations), but due to 
them first calculating static spin susceptibility in momen- 
tum space 

y(q) - V 7W , [E^^jp + q)] ~ np [i^.(p)(U 



(we shouldn't worry here what the matrix element AA is) 
and then calculating x (Ry ) making a Fourier transfor- 
mation 



X(R..)= / d\x{ci)e^'^-^^^. 



(50) 



Both integrals turn out to be ultra-violet divergent, and 
cut-offs should be introduced. We, on the other hand, 
calculated directly x in real space representation, thus 
avoiding the problem of divergence of the integrals com- 
pletely. 

There is another problem with calculating the RKKY 
interaction (in normal metals) which has a long history 
2^ . 23| ; it arises when we combine the integrals and 
([50|) into a single double integrals. The problem is which 
integration: with respect to q or with respect to p we 
should do first. We also avoid this problem completely. 

The contact exchange interaction we used can be eas- 
ily justified in the case of s-wave orbital of the magnetic 
adatom [llj] . The case of d-wave orbitals is more compli- 
cated. To find the physically meaningful form of Kondo 
perturbation, it is appropriate to go back to the possi- 
ble origin of the Kondo model, i.e., the Anderson model. 
Following seminal paper by Schrieffer 2j|, let us spec- 



ify the magnetic impurity as being the S-state ion, say 
Mn++, whose d-shell has the configuration 5*5/2. Since 
the S'-state ion cannot change the orbital angular mo- 
mentum of a conduction electron, one should use states 
which transform according to the irreducible representa- 
tions of the point group of the crystal about the impurity 
center [24]. 

Such approach for the case of d-wave orbitals was re- 
alized by Zhu et al. [1^ (see also Ref. [HI). Considering 
the magnetic impurity above the center of the honeycomb 
(plaquette impurity), they started from the classification 
of the degenerate 3d-orbitals of the magnetic atom with 
respect to irreducible representations of the symmetry 
group Cgt, and inferred that d^^ belongs to Ai, {dxz,dyz) 
belong to Ei and (dx'i -y^ ^ dxy) belong to E2 representa- 
tions. Specifying their approach, we'll take into account 
hybridization of the d-orbitals of the magnetic impurity 
with the states of the carbon atoms around the pla- 
quette. The selection rules for matrix elements demand 
that from the states \i >, where i G V, and V is the set of 
sites surrounding the plaquette, we'll chose combinations 
realizing the same representations as above. Thus the 
hybridization Hamiltonian for the 3d magnetic impurity 
in terms of the irreducible reps of the system will take 
the form 



(51) 



X,a,ie'P 



where operators (/) create (annihilate) electrons at 
the d-orbitals of the magnetic impurity, and index A enu- 

From Eq. 



merates the orbitals dz^ , dxz , dyz , d^^ 



-y^ 7 "-xy 
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(I5T|) . following Ref. 2j| under appropriate assumptions 



we can get the p — d exchange model 
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Hpd — — 



X,a,l3,i,j£'P 



c 



(52) 



In Ref. m , the p — d exchange Hamiltonian (for the 
so-called coherent case) was previously taken in a very 
specific form 

Hpd = -J ^ S-ct„^4„Cj73, (53) 

which in fact takes into account only the hybridization 
between dz'^ and the combination of the p-states on the 
plaquette, realizing irreducible representation Ai, that is 



E 



> /\/&. Such specific form led to the conclu- 
sion that 1/|R — R'l'^ term in the RKKY interaction be- 
tween the plaquette impurities vanishes. When the gen- 
eral form of the hybridization Hamiltonian ()5ip is taken 
into account, this conclusion seems to us unjustified. 
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